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ABSTRACT 

We use 1 kpc resolution cosmological AMR simulations of a Virgo-like galaxy cluster 
to investigate the effect of feedback from supermassive black holes (SMBH) on the 
mass distribution of dark matter, gas and stars. We compared three different mod- 
els: (i) a standard galaxy formation model featuring gas cooling, star formation and 
supernovae feedback, (ii) a "quenching" model for which star formation is artificially 
su ppressed in massive halo s and finally (iii) the recently proposed AGN feedback model 
of iBooth fc Schavi (|2009D . Without AGN feedback (even in the quenching case), our 
simulated cluster suffers from a strong overcooling problem, with a stellar mass frac- 
tion significantly above observed values in M87. The baryon distribution is highly 
concentrated, resulting in a strong adiabatic contraction (AC) of dark matter. With 
AGN feedback, on the contrary, the stellar mass in the bright central galaxy (BCG) 
lies below observational estimates and the overcooling problem disappears. The stellar 
mass of the BCG is seen to increase with increasing mass resolution, suggesting that 
our stellar masses converges to the correct value from below. The gas and total mass 
distributions are in better agreement with observations. We also find a slight deficit 
(~10%) of baryons at the virial radius, due to the combined effect of AGN-driven con- 
vective motions in the inner parts and shock waves in the outer regions, pushing gas 
to Mpc scales and beyond. This baryon deficit results in a slight adiabatic expansion 
of the dark matter distribution, that can be explained quantitatively by AC theory. 

Key words: black hole physics - cosmology: theory, large-scale structure of Universe 
- galaxies: formation, clusters: general - methods: numerical 



1 INTRODUCTION 

Galaxy clusters are ideal laboratories to study galaxy for- 
mation in a dense environment. Galaxies in clusters are ob- 
served in many morphological types, from blue extended spi- 
rals to red massive elliptical spheroids. The origin of the 
morphological evolution of galaxies in clusters is still poorly 
understood. Tidal and ram pressure stripping trigger fast 
evolution in the properties of accreted satellites, while com- 
plex gas cooling and heating processes control the amount 
of stripped gas that can actually reach the central region 
of the cluster. In this context, the origin of bright cluster 
galaxies (BCG) still raises many questions. In the current 
cosmological framework, the formation of galaxies at the 
bright end of the luminosity function is still affected by the 
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so-called "overcooling problem": using both semi-anaytical 
models and computer simulations, it has been shown that 
the massive galaxies are predicted to be too bright and too 
blue when compared to massi ve galaxies in the nearby un i- 
verse (see the recent review bv lBorgani fc Kravtsovl l|2009l )'). 
As a consequence, these models find a stellar content in mas- 
sive clusters that is s ignificantly above the observed values 
jKravtsov et al.|[2005h , even including rathe r extreme super- 
novae feedback recipe (|Borgani et al.|[2004i '). 

In order to overcome this issue, feedback from su- 
permassive black holes (SMBH) have been proposed as a 
mechanism to prevent gas from accumulating in the clus- 
ter core. The so-called AGN feedback scen ario has received 
supp o rt from theoretical con s iderations (iTabor fc BinnevI 
1 19931 : ICiotti fc OstrikeJ 1 19971 : ISilk fc Reei 1 19981 ). but "the 
strongest evidence comes from the correlated observations 
of X-ray cavities and radio blobs in massive clusters. 
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These features are usually interpreted as buoyantly ris- 
ing bubbles of high entropy material injected in the clus- 
ter core by jets of relativistic partic le s. Detailed mod- 
els of bubbles JChurazqv et al] I2001I: Ruszkowski et al.) 
20041: iBriiggen et al l I2OO5I) and iets llRevnolds et al.ll200ll : 



Omma et all l2004l : I Cattaneo fc Tevssieil 120071 7 have been 
proposed in the context of cluster cores heating, usually 
based on spherically or planar symmetric, idealized cluster 
configurations. Based on simple energetic arguments, it is 
possible to relate the growth of SMBHs at the centre of mas- 
sive galaxy sphe roids to the energ y required to unbind the 
overcooling gas l|Silk k, Reeslll99g ). The idea of star forma- 
tion being regulated by AGN feedback at the high mass end 
of the galaxy mass function has been applied first quite suc- 
cessfuUy to hydrody namical simulations of galaxy merger 
IIMatteo et al. I I2OO5I ) and then to semi-analytical models 
of galaxy formation llCroton et al. 2006: Bower ct al. 2006.: 



Cattaneo et~ai]|2006l : lLucia fc Blaizotll2007i : ICattaneo et alj 
20091 ). 



AGN feedback models have been included only recently 
in co smologic al simulations of galaxy groups and clusters 
llsiiacki et a l.l l2007l : iPuchwein et al.1 l2008l : iMcCarthv et al.1 
I2OO9I ). Although the d etailed physical modehng of SMBHs 
grow th usually diflFer l|Siiacki et al.1 l2007l : iBooth fc Schavd 
I2OO9I ). these simulations, exclusively based on the GAD- 
GET code (Springell lioosh . have obtained very encourag- 
ing results re garding the global p r operties of the simu- 
lated clusters ifpuchwein et al.ll2008l: iMcCarthv et al.ll2009l : 
iPuchwein et al.l I2OI0I : iFabian et ail boid ). In this paper, 
we report on the first AMR high-resolution simulation of 
a Virgo-like cluster, following both SMBH and star for- 
mation, with the associated feedback. A companion paper 
is expl oring the propertie s of a jet-based, AGN feedback 
model (iDubois et al.(l2010l '). We use the AMR code RAM- 
SES, which differs significantly from the other code used so 
far to model AGN feedback in a fully cosmological simula- 
tion, namely the GADGET code. Although AMR schemes 
sufi^er from larger advection errors than SPH, which is a 
strictly Galilean invariant method, there are some definitive 
advantages of using AMR in this context: although the total 
energy (kinetic -I- thermal -I- potential) is conserved only at 
the percent level, it is strictly conserving for the fluid energy 
(kinetic -|- thermal), which is very important in presence of 
strong shocks, and it captures h ydrodynamical instabilities 
more realisticallv llAgertz et al.1 12007: Wadsley et al., ,2008,: 
iMitchefl etailbOOg T which is very important in presence 
of convective motions of buoyant gas. It also captures gas 
stripping of infalling satellites due to hydrodynamical insta- 
bilitie s more realistically than standard implementations of 
SPH (|Agertz et al.ll2007l ). 

We have adapte d the SMBH growth model of 
iBooth fc Scha ve' (2009') to the sink particle method for AMR 
presented by Krumholz ct al. (2004)^_With respect to the 
previous work of Boot h fc Schavd ((20091) and follow-up pa- 
pers, we have made significant improvements over the OWL 
simulations suites in term of mass and spatial resolution, but 
only for one sin gle zoom-in simu latio n of a Virgo-like cluster . 
Recentlv. iPuchwein et al.1 l|2008t ) and IPuchwein et al.1 (I2OI0I ) 
have also used the GADGET code, but a different AGN 
feedback model, to perform zoom-in simulations of groups 
and clusters of galaxies, with however a mass resolution and 
a gravitational softening length not as good as the one we 



used in our high resolution case. Note also that the minimum 
smoothing length in these SPH simulations is usually much 
smaller than the gravitational softening length. In this pa- 
per, we have specifically chosen a Virgo-like cluster, in order 
to compare our results with the very detailed observations 
that are available for this well-known astronomical object. 
We will focus here on the mass distribution of the three main 
components, namely dark matter, gas and stars. 

The paper is organized as follows: the first section is 
dedicated to numerical methods and physical ingredient 
(cooling, star formation and AGN feedback), while our sec- 
ond section presents our results, comparing our three mod- 
els. The final section is left for discussion. 



2 NUMERICAL METHODS 

In this section, we describe the numerical techniques and the 
initial conditions we used to model our Virgo-like cluster. As 
it is now customary for cosmological simulations, we used a 
zoom-in (or volume renormalisation) technique to focus our 
computational resources on a specific region of a 100 Mpc/h 
periodic box. We adopt a standard ACDM cosmology with 
= 0.3, = 0.7, »i, = 0.045 and i/p = 70 km/s/Mpc. 
We have adopted the lEisenstein fc Hul lll998[) tra nsfer func- 
tion and the graf ic package ('BertschinEer"200l') in its par- 
allel implementation mpgraf ic (Prunot ct al. 2008) to gen- 
erate our initial conditions. From a first low resolution dark 
matter only run, we built a catalog of candidate halos whose 
virial masses lie in the range 10" to 2 X 10" Me/h. From 
this catalog, we have extracted our final halo based on its 
mass assembly history: its final mass 

is already in place around z = 1, so it can be consid- 
ered as a well relaxed cluster hy z — 0. The final halo 
mass has been measured to be M200C = 1.04 x lO'** Mq 
or M500C = 7.80 X lO''' Mq, where indice c refers to the 
critical density. 



2.1 Simulation parameters 

We have performed two simulations, one at low resolution, 
for which the initial grid effective size was 1024'^, and one 
at high resolution, with effective grid size of 2048^. From 
this initial grid, we have extracted high resolution particles 
only in the Lagrangian volume of the halo, and we have used 
larger and larger mass particles to sample the cosmological 
volume, so that the total number of particles in the box was 

5.2 X lO'' (resp. 22 x 10*^) for only 2.6 x 10^ (resp. 19 x 10*^) 
in the central region, and 10** (resp. 8 x 10**) inside the fi- 
nal virial radius for the low resolution (resp. high resolu- 
tion) simulation. The dark matter particle mass is therefore 
6.5 X 10'' (resp. 8.2 x 10^) Mq/Ii and the baryons resolution 
element mass is 1.2 x lO'' (resp. 1.4 x 10^) Mo/h. 

The AMR grid was initially refined to the same level 
of refinement than the particle grid (1024^ and 2048'^), and 
7 more levels of refinements were considered. We imposed 
the spatial resolution to remain roughly constant in phys- 
ical units, so that the minimum grid cell stayed close to 
Axmin = L/2*""'- with ^max =17 (rcsp. 18) at 2 = 0. We 
therefore reached a spatial resolution of Axmin — 1 kpc for 
the low resolution simulation and Axmin — 500 pc for the 
high resolution one. The grid was dynamically refined up to 
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Run 




^gas 


nit 








in 10*^ Mq 




in kpc/h 


Low res 


65 


12 


2.4 


0.76 


High res 


8.2 


1.4 


0.3 


0.38 



Table 1. Mass resolution for dark matter particles, gas cells and 
star particles, and spatial resolution (in physical units) for our 2 
sets of simulations. 



this resolution using a quasi-Lagrangian strategy: when the 
dark matter or baryons mass in a cell reaches 8 times the 
initial mass resolution, it is split into 8 children cells. We 
reached z = with 14 x lO'' (resp. 68 x lO'') cells for the low 
(resp. high) resolution run, including split cells. 

Gas dynamics is modeled using a second-order un- 
split Godunov schem e (|Tevssieill2002[ : iTevssier et al.ll2006l : 
[Fromang ct al.' '2006*) based on the HLLC Riemann solver 
( Toro e t al. 1994). We assume a perfect gas equation of state 
(EOS) with 7 — 5/3. Gas metallicity is advected as a pas- 
sive scalar, and is self-consistently accounted for in the cool- 
ing function. We al so considered the s t andar d homogeneous 
UV background of'Haardt & Madad !\199di ). but we modi- 
fied the starting redshift, extrapolating the average intensity 
from Zrcion = 6 to Zroion ~ 12, in order to account for the 
early reionization expected in such a proto-cluster regions 
l|lliev et al.ll2008l ). 



2.2 Galaxy formation physics 

As gas cools down and settles into centrifugally supported 
discs, we need to provide a realistic model for the interstellar 
medium (ISM). Since the ISM is inherently multiphase and 
highly turbulent, it is beyond the scope of present-day cos- 
mological simulations to try to simulate it self-consistently. 
It is customary to rely on subgrid models, providing an ef- 
fective EOS that capture the basic turbulent and thermal 
properties of this gas. Models with various deg rees of com- 
plexi t y have been proposed in the literature ("Yepcs ct al. 
Il997l : ISpringel fc Hernauistil2003l : ISchave fc Vecchia ,2008.) . 
We follow the simple approach based on a temperature floor 
given by a polytropic EOS for gas 



Tfloor = T* I — 5- 



(1) 



where n* =0.1 H/cc is the density threshold that defines 
the star forming gas, T* = 10* K is a typical temperature 
mimicking both thermal and turbulent motions in the ISM 
and r = 5/3 is the polytropic index controlling the stiffness 
of the EOS. Gas is able to heat above this floor, but cannot 
cool down below it. Note that because of this temperature 
floor, the Jeans length in our galactic disc is always resolved. 
We also consider star formation using a similar phenomeno- 
logical approach. In each cell with gas density larger than 
n*, we spawn new star particles at a rate given by 



■ e* - 



with ts 



37r 
32Gp 



(2) 



where tg is the free-fall time of the gaseous component and 
e, = 0.01 is the star formation efficiency. The star particle 



mass depends on the resolution and was chosen to be 2.4 x 
10^ (resp. 3x 10^) M© for the low (resp. high) resolution run. 
For each star particle, we assume that 10% of its mass will 
go supernova after 10 Myr. We consider a supernova energy 
of 10^^ erg and one Mq of ejected metals per 10 Mq average 
progenitor mass. This supernovae feedback was implemented 
in the RAMSES co de using the "delayed cooling" scheme 
(|Stinson et al.l l2006Y 

Up to this point, we use rather standard galaxy for- 
mation recipe, which have proven only recently to be 
quite s uccessful in reproducing the properties o f field 
spirals (iMayer et all l2008l : iGovernato et al.l |2009| . l2010l : 
lAgertz et al.ll2010l ). These recipe have been shown to re- 
produce basic galaxy properties like Kennicutt-Schmidt 
law, s tar formation rates, ga l actic winds foubois fc Tevssieil 
I2OO8I : iDevriendt et al.ll2010l : lAgertz et al.| |2010l) . The same 
recipe are only marginally su ccessful when one consid- 
ers small grou ps (iFeldmann et al.ll2010l), but they fail on 
cluster scales |Borgani et al.l 20041 : iKravtsov et al.l l2005l : 
iBorgani fc Kravtsovl 20091 '). In the present paper, in order to 
check that our results are compatible with previous work, 
and to set a reference point, we have performed simulations 
with only standard galaxy formation physics (labelled SF 
runs in the foUowings). 



2.3 Star formation quenching 

The main problem we have to face in standard cosmolog- 
ical simulations on cluster scales is the striking excess of 
mass locked into stars. Using rather extreme stellar feedback 
models, previous authors report that the simulated stellar 
mass fraction lies in the range 35 - 60% , depending on the 
cluster mass (|Borgani fc Kravtsovll2009l ). Since only 10% of 
the baryonic mass is observed in the galaxies, this would 
require a large amount of stars hidden in a diffuse compo- 
nent such as the intracluster light (ICL). This last scenario 
is however severely co nstrained by observations of the ICL 
(|Gonzalez et al.|[2007l ') and does not appear to be plausible. 

There is growing theoretical and observational evidence 
that star formatio n is shutdown above a critical halo mass 
Mc ~ 6x 10" MM l|Cattaneo et al.lliooel 'l. lBirnboim fc Dekell 
((2003i) have proposed that this critical mass is related to the 
stability of accretion shocks, and to a transition from cold 
to hot mode of gas accretion. Although this transition in 
the nature of the acc retion flows has been confirmed by nu- 
merical simulatio ns (iKeres et al ] I2OO5I : lOcvirk et al. I l2008l : 
iDekel et al.ll2009l ). the simulated star formation was not 
observed to decrea se significantly above the critical mass 
(|Ocvirk et al.ll2008l ). This might be due to insufficient mass 
and spatial resolution, so that heating processe s, not prop- 
erly c aptured, cannot balanc e radiative cooli n g ([Naab et al.l 
l2007h . On a different side, ICattaneo et all (120061 ) argued 
that this modification of the halo properties may create fa- 
vorable conditions for AGN feedback to be effective, but 
AGN feedback is still needed t o prevent the overcoolin g 
problem above the critical mass (jPekel fc Birnboimll200d '). 
Without specifying the underlying heating process, the crit- 
ical mass argument boils down to stopping (or quenching) 
gas cooling and star formation in the central galaxy for halo 
masses larger than ~ lO'^^ Mq. 

In this paper, in order to test this hypothesis, we have 
used a simple phenomenological model to quench star forma- 
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tion in massive galaxies. Since our standard model (SF runs) 
obviously suffers from a strong overcooling problem, we need 
to actively suppress gas cooling and star formation above 
the critical mass. If the stellar mass density is greater than 
0.1 H/cc, and if the stellar 3D velocity dispersion is greater 
than 100 km/s, we switch off star formation and gas cool- 
ing. In this way, star formation in discs is unaffected, since 
the velocity dispersion is smaller than the chosen threshold. 
On the other hand, star formation is suppressed in massive 
spheroids. The main advantage of this quenching model is 
its simplicity: a l thoug h it captures the scenario proposed by 
ICattaneo et al.l (|2006l ) . it does not require any complex AGN 
fee dback mode l , nor t he mass and spatial resolution reached 
bv lNaab et al.l (20o3) in their early-type galaxy simulation. 
However, as we demonstrate in this paper, this simple ap- 
proach doesn't solve the overcooling problem in our Virgo 
cluster simulation. 



2.4 SMBH growth and associated feedback 

Beside our standard galaxy formation and our quenching 
scenario simulations, we explore a model for which SMBH 
growth and feedback is considered. In a nut shell, our recipe 
is based on the sink particle met hod for grid-based codes 
designed bv lKrumholz et all (120041'). with th e AGN feedback 
model proposed bv iBooth fc Schavd ()2009h . We shall now 
briefly summarize these techniques. 



2.4-1 Seed particles 

In the two main competing SMBH formation mo dels: slow 
growth from Pop HI stars (jMadau fc Reed 200ll) or direct 



collapse of a low angu lar momentum halo ijBromm fc Loebl 
I2OO3I : iBegelman et al.l | 2006l ). the seed SMBH is believed to 
grow quickly to 10^ Mq, before starting to affect its environ- 
ment and enter the self-regulated regime. Although the ques- 
tion of intermediate mass black holes is still vigorously de- 
bated, this characteris tic ma ss is often considered as the ini- 
tial seed SMBH mass (|Li et a l. 2006; Pclupessy ct al. 2007), 
since it is one order of magnitude smaller t han the smallest 
SMBH observed on the M bh — cr relation ijGebhardt et al.l 
I2OOOI : iGiiltekin et ahllioogl ). At each main time step during 
the course of the simulation, we identify potential candidate 
regions for seed black hole formation using the following cri- 
teria: 

• the stellar density has to be greater than 0.1 H/cc 

• the stellar 3D velocity dispersion has to be greater than 
100 km/s 

• the gas density has to be greater than 1 H/cc 

• no other sink particle is present within 10 kpc 

If these four conditions are fulfilled, we create a sink particle 
of fixed mass Afan ~ 10^ M© . These particles represent our 
seed black holes. Note that our third condition ensures that 
seed black holes form in the nuclear region of star forming 
disks, where the gas density is probably much larger than 
1 H/cc. Because our limited resolution, we cannot choose an 
arbitrary high density threshold, otherwise SMBH will never 
form. On the other hand, because star formation is a very in- 
efficient process, large enough galaxies devoid of SMBH will 
always reach this gas density threshold and trigger SMBH 



seeding. Our second condition requires a minimum line- 
of-sight velocity dispersion o\d — 6 km/s, in agreement 
with the observ ed Mbh — o relation (|Gebhardt et al.ll2000l : 
IGiiltekin et al.l 2009). Our last condition ensures that no 
new seed SMBH will be created in a galaxy that is already 
hosting an old SMBH (or at least within 10 kpc from its 
center) . 



2.4.2 Svnk particle evolution 

Each black hole i s cons idered as a sink particle, as defined in 
iKrumholz et al.l (|2004l ). We recall briefiy the method here. 
We consider around each black hole a sphere of fixed physical 
radius rsink = 4 Ax, where Ax is our spatial resolution in 
physical units, so that r^ink — 4 kpc (resp. 2 kpc) for the 
low (resp. high) resolution run. We assume that the SMBH 
mass distribution inside the sink is homogeneous, and this 
homogeneous sphere is added to the total mass density when 
solving for the Poisson equation. The sink particle is then 
advanced in time by interpolating the gravitational force 
back to the sink position using the inverse CIC scheme. For 
each sink, we compute the Bondi-Hoyle accretion rate 



Mb 



ttboost 



47rG^MgH/o 

(c2 +tt2)3/2 



(3) 



where p, Cs and u are the average gas density, sound speed 
and relative velocity within the sink radius. These aver- 
age CQWiitittes_ax^ following the scheme proposed 
bv IKrumholz et al] (12004 ). The parameter Qboost was in- 
troduced by Springel et al.l (|2005l ) to account for unresolved 
multiphase turbulence in the SMBH environment. Although 
this parameter was fi rst chosen constant at aboost — 100 
l|Springel et al.ll2005h . iBooth fc Scha^ Hooi) argued that 
this parameter should be close to one in low density regions, 
while its value should increase in high density regions, in 
order to match the subgrid model used for the unresolved 
turbulence in the disks. Based on extensive numerical experi- 
ments, they proposed the following phenomenological model 
to boost Bondi-Hoyle accretion. 



Q^boost 
Q^boost 



n* 

1 otherwise 



if riH > n« = 0.1 H/cc, 



(4) 



Note that this model has no real physical justification, and 
that it depends crucially on the underlying EOS model. 
We have been very careful in using this boost factor in 
conjunction with the s ame EOS model (see Equ. [T]) as in 
iBooth fc Schavd (|200S " 



1). 

As advocated bv ISpringel et al.l (|2005l ). the accretion 
rate onto the SMBH cannot exceed its Eddington limit given 
by 



Med = 



47rGMBHmp 



with 



0.1, 



(5) 



so that the final accretion rate is computed using Afacc ~ 
min(AfBH, A^ed). At each time step, a total gas mass of 
Aifacc At is removed from all cells within the sink radius, with 
the same weigh ting scheme as the on e used to define aver- 
age quantities l|Krumholz et al.ll2004l ). In order to prevent 
the gas density to vanish or become negative, we allow a 
maximum of 50% gas removal at each time step. 
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Figure 1. Time evolution of the accretion rate of the largest 
SMBH in the cluster. The Eddington limit is shown as the red 
line, the corresponding SMBH mass indicated on the right axis. 
The high resolution simulation finds a similar (within a factor 
of 2) SMBH mass (shown in blue). We clearly see two modes of 
accretion, with episodic bursts for which the SMBH accretes at 
or close to the Eddington limit, and long periods in between for 
which the SMBH accretion rate is at or close to the Bondi rate. 



24.3 AGN feedback 

In the proposed model for SMBH growth, a key ingre- 
dient is the associate d feedback mo d el. As demonstrated 
by many auth ors (Sii acki et al.ll2007l : ICattaneo fc TevssieJ 
I2OO7I : iBooth fc Schavei l2009h . the Bondi-Hoyle accretion 
model allows the SMBH growth history to be self-regulated 
by injecting thermal energy in the surrounding gas: if the 
SMBH mass is too small, the associated feedback will be in- 
efficient and the surrounding gas will remain cold and dense, 
boosting the accretion rate up to the Eddington limit. The 
SMBH mass will grow exponentially fast, with e- folding time 
scale equal to the Salpeter time ts — 45 Myr. As soon as the 
black hole mass is large enough, feedback processes heat and 
eventually unbind the surrounding gas, so that the accretion 
rate drops and becomes Bondi-Hoyle limited. This bimodal 
behavior is illustrated in Figure[Tl where one can see the time 
evolution of the accretion rate of the most massive SMBH in 
the simulated region. Short bursts of Eddington limited ac- 
cretion are followed by long quiescent epochs of Bondi-Hoyle 
limited accretion, self-regulated by SMBH feedback. 

In order to allow for this self-regulated SMBH growth, 
efficient feedback schemes are mandatory. Although we 
do see clear sign atures of AGN feedback in clus ters 
jArriaud et al.l [l98 4: Fabi an et al.l |2000; McNamara ct~al] 
|2001^ ■ the physical processes at the origin of this 
energy injection are still unclear: radiative feedback 
JCiotti fc Ostrikeij|200ll'). cosmic rays iBrtiggen et al.ll2002l : 
ICha ndran fc Raseral 120071 ') or strong shocks (see the review 
of iBegolman C200j)). A common property of these various 



models is that they require a very good spatial resolution 
to be captured realistically in hydrodynamical simulations. 
The most advanc ed modeling so f ar ha s been using AGN- 
driven bubbles jChurazov et al.l I2OOII : iRuszkowski et ahl 



Qmma et al. 2004; Cattaneo fc Tevssiej 20071 ). leading to 
turbulent co nvective mo tions and "shockle ts" escaping the 
clust e r core ( Chandra nfc Rasera„2007. : .Rasera fc ChandranI 



I2OO8I : ISharma et al.l |2009| ). In some case, depending on 
the injected energy, these AGN-driven flows can drive 
strong shock wav es that can travel to very large distances 
(Baldi et al."2009'). In current cosmological simulations, nu- 
merical resolution does not allow these effects to be self- 
consistently modeled. As it is customary under these cir- 
cumstances, we rely on a more phenomenological model. 

In the cosmologic al context, the model proposed by 
iBooth fc Sciia^ l|2009D a ppears to be easier to implement 
than the one proposed bv lSiiacki et al.l (|2007l ). Moreover, it 
relies on only one main free parameter, the coupling effi- 
ciency Ec, that can be calibrated to the observe d Mbh — o" 
relati on to the fiducial value ec — 0.15 (Boot h fc Schavj 
I2OO9I ). We have adapted their model to the RAMSES code, 
using the following approach: at each time step, we com- 
pute the SMBH feedback energy as a fixed fraction of the 
rest mass energy of the accreted gas, multiplied by the "cou- 
pling efficiency" e^, 



AE = EcErA'/a. 



'At. 



(6) 



This energy is not released immediately in the surrounding 
gas. It is instead accumulated over many time steps and 
stored into a new SMBH related variable Eagn, so that we 
can avoid atomic line cooling to radiate th is energy instan- 
tanoously. Following the trick proposed bv lBooth fc Schavei 
(2009), we release this accumulated energy inside the sink 
radius when 



-Eagn > g'^'s 



(7) 



where mgas is the gas mass within the sink radius and Tmin 
is the minimum feedback temperature. As soon as Tmin is 
chosen above lO'^ K, the critical temperature below which 
metal line cooling becomes very efficient, the resulting feed- 
back scheme does not depend on the chosen value for Tmin. 
We adopt here Tmin ~ K. As can be seen on the previous 
equation, this threshold energy depends directly on the gas 
density in the environment of the black hole. When dense 
and cold gas is present, more energy is required to reach the 
threshold. After enough mass has been accreted, a strong 
burst of energy is released, that will unbind the surrounding 
dense gas. On the other hand, when only diffuse, hot gas is 
present, the threshold energy is much easier to reach, and 
feedback proceeds in a quasi-continuous fashion. In some 
sense, based on this rather simple recipe, we can account for 
both the "quasar mode" and the "rad io mode" of the AGN 
feedback model of ISiiacki et al.1 (|2007l ). 



3 RESULTS 

In this section, we describe the properties of our simulated 
cluster for the 3 different models, labelled "SF", "quenching" 
and "AGN" in most of the Figures. We will focus our analy- 
sis in the final mass distribution, and compare, whenever it 
is possible, to actual Virgo cluster data. We present mostly 
low resolution data, although we also compare low and high 
resolution results to discuss convergence properties. 
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Figure 2. Maps of the mass weighted temperature of the simulated cluster at four diflfcrent redshifts. Units are in keV. The first redshift 
in the series corresponds to the strong SMBH outburst seen in Figure [T] around a ~ 0.62. 



3.1 SMBH Growth and Feedback 

In Figure [TJ we show the accretion rate of the most massive 
SMBH as a function of expansion factor. This plot is rele- 
vant only for the "AGN" simulation. Also shown is the Ed- 
dington accretion rate, directly proportional to the SMBH 
mass. It appears quite clearly in this plot that the most mas- 
sive SMBH grows discontinuously, during very short Edding- 
ton limited accretion events, or, at late time, by accreting 
other black hole s, in good agreemen t with semi-analytical 
predictions from lMalbon et al] (|2007l ). In the low resolution 
simulation, the final mass reaches Mbh ~ 2.1 x 10^ M© 
after a last merger around a ~ 0.85. In the high resolu- 
tion simulation, the SMBH mass is twice as large, with 



a slightly different evolution, and reaches the final value 
of AfsH ~ 4.2 X 10^ M0. The SMBH mass in M87 has 
been est imated around 4 x 10 ^ M^ using dynamical con- 
straints (iMacchetto et"allll997l : iGebhardt fc Thomas! [2OO9I : 
iGiiltekin et al.ll2009l l. in very good agreement with our high 
resolution predict ion. Note that M87 SMBH is close to the 
Mbh-o" relation l|Giiltekin et al. I [2009. 1. Since the free pa- 
rameter tc was calibrated to tc = 0.15 bv iBooth fc Schava 
(j2009l ) on the Mbh-o" relation, our AMR simulation appears 
to be consistent with their SPH results. 

The SMBH activity, quite strong before z = 1, declines 
slightly until the present epoch. In our simulation, this early 
activity is due to an early phase of frequent mergers feeding 
the SMBHs very efficiently. Strong and repeated outbursts of 



AGN and Mass Distribution 7 



199ll : 
20081 ). 



energy are launching strong shock waves in the IGM, rising 
the IGM entropy within the whole proto-cluster region. This 
early epoch can therefore be considered as representative of 
the pre-heating scenario advocated by several authors to ex- 
plain structural properties of galaxy clusters |Kaiserl 
iPonman et al] 1 19991 : iBabul et al.l |2002| : iDave et al.l \2 
On the other hand, at later epochs (2 < 1), when the clus- 
ter mass is finally assembled, AGN feedback prevent gas 
from overcooling and from accumulating in the core. This is 
well illustrated by the sequence of temperature maps shown 
in Figure [2l just after the strong AGN outburst occurring 
at a ~ 0.62 (see Fig. [T]). The first image show the mass- 
weighted projected temperature within the whole cluster, 
just at the time of the outburst. Slightly after (in the sec- 
ond frame), the whole cluster has been significantly heated, 
with buoyantly driven plumes of hot gas escaping the cluster 
core JChandran fc Rascra 2007; Cattanco & Tcyssicr 2 0OT| : 
iRasera fc Chandranll2008l : ISharma et al.ll2009l ). In the next 
frame, a strong shock, visible as a sharp temperature dis- 
continuity, develops close to and beyond the virial radius. 
The last frame shows the cluster back to hydrostatic equi- 
librium, waiting for the next AGN outburst. Note that in 
the core region, where the density is high enough for X-ray 
detections, only buoyantly rising bubbles are seen. In our 
simulation, shock waves form only in the outer part of the 
cluster, with a Mach number of a few. Because of the rather 
low gas density, we believe that their detection is quite chal- 
lenging, e xplaining why there is so little evidence of their 
presence jBourdin et al.ll2004 ). 



Although the feedback recipe doesn't explicitly account 
for it, the AGN energy deposition typically proceeds in 
two modes: a strong, energetic mode or "quasar mode", 
which in our case corresponds to the Eddington luminos- 
ity and reaches 5 x 10^^ erg s~^, and a quiescent mode 
or "radio mode", with a Bondi-Hoyle-limited luminosity of 
5 X 10*^ erg s~^. If we now take into account the coupling 
efficiency parameter ec — 0.15 in our energy estimate, using, 
from Figure [1] Mace — 10~^ Mq/jt into Equation [S] we ob- 
tain a total luminosity of 9 x lO*'' erg s~^ during radio mode. 
Although the X-ray luminosity of the active nucleus of M87 
observed with Chan dra is I/x, 0.5-7 kcv ~ 7 x 10*'' erg s~^ 
l|Matteo et al.|[2003l l. quite close to our total luminosity in 
the radio mode, the mechanical power of the observed jet 
is sig nificantly larger Pjet — 3.3 x 10*'^ erg s~^ (|Allen et al.l 
l2006i V This is consistent with observational estimates of the 
Bondi accretion rate in M87, around Mace — 10~^ M0/yr 



l|Matteo et al.ll2003l: lAUen et al.ll2006l 'l and a coupling effi- 
ciency Ec ^ 0.2 (|Allen et al.ll2006l '). This suggests that M87 
SMBH lies in an intermediate state, between our radio and 
quasar modes. Nevertheless, the fact that our final black 
hole mass match well observational constraints of M87 is of 
great importance: it means that the available rest mass en- 
ergy that the black hole can release into the forming cluster 
over its entire history is consistent with M87. This doesn't 
mean that the actual amount of energy deposited into the 
X-ray emitting gas today is correct. With these limitations 
in mind, we can now predict the global properties of our 
simulated Virgo cluster, in particular the mass distribution 
in stars, gas and dark matter. 




100 
r(kpc) 



1000 



Figure 4. Cumulative stellar mass profile for our simulated clus- 
ter at z = 0. Tlie blue line i s tlie stellar mass profile of M87 from 
iGebhardt &: TliomasI | |2009|) . 
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Figure 5. Gas density profile for our simulated cluster at 2 = 0. 
The blue line is the fit of the gas dens i ty usi ng the X-ray emissivity 
profile of M87 from lChurazov et al.l | |2008|'I . 



3.2 The mass distribution of stars and gas 

We have plotted in Figure |3] the surface brightness of our 
simulated Virgo cluster in the SDSS i band. One clearly 
sees many satellite galaxies orbiting around the BCG and 
a rich structure in the intracluster light component. In the 
SF case, the BCG appears as a very bright, gas rich and 
disc-like object. We have also plotted in Figure|4]the stellar 
mass profile from this SF simulation: we see immediately 
that the BCG stellar mass (measured at 20 kpc from the 
centre) is a factor of 10 to large, when compared to the 
observational estimate of lCebhardt fc Thomaj ()2009( ). This 
is the classical result of the over-cooling problem that oc- 
curs for cluster simulations with standard galaxy formation 
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Figure 3. Maps of the projected stellar luminosity (left panels, units are I band absolute magnitude) and gas mass weighted density 
(right panels, units are in H per cc) in the simulated cluster at z = for our three models. 
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Run 


ji rtot 
*-'200c 




/"bar 


/. 


/gas 


SF 


1.2 X lO^-* 


Mq 


16% 


8% 


8% 


Quenching 


1.2 X 10^-* 


M0 


16% 


3% 


13% 


AGN 


1.1 X 10" 


Mq 


13% 


1% 


12% 



Table 2. Mass fractions inside the virial radius for our three 
different models. The universal baryon fraction we used in this 
paper is 15%. 



physics. Our second scenario, the Quenching run, partially 
alleviates this problem. As can be seen on the surface bright- 
ness maps (Fig. [Sjl , the BCG and the most massive satellites 
are now much dimmer. The total stellar mass within 100 kpc 
is in much better agreement with observations, although still 
slightly larger. The stellar mass profile, shown in Figure IH 
is however significantly different. When one looks now at 
the gas distribution in the cluster, it becomes quite obvi- 
ous that this Quenching scenario is far from being a viable 
solution. We have plotted in Figure O the mass- weighted, 
projected gas density. We see a massive gas clump in the 
cluster core for both the SF and the Quenching runs. For 
sake of comparison, we have plotted the simulated gas den- 
sity profil es for our simula t ions a nd the best-fit /3-modeI for 
M87 from lChurazov et al.1 l|2008l ). The gas density is the SF 
run is a factor of 100 too large in the core of the cluster, 
and it is even worse for the Quenching scenario, by an ad- 
ditional factor of 2. We therefore conclude that for both SF 
and Quenching models, the simulated cluster suffers from a 
strong overcooling problem, with the build up of a dense, 
concentrated BCG, for which the stellar mass or the gas 
mass (or both) are in far in excess of those observed in M87. 

We now turn to the analysis of our AGN model. The 
stellar surface brightness map is by far the dimmest of our 
3 models: star formation has been dramatically reduced, 
even more than our Quenching scenario. The stellar mass 
profile is now below the observational constraints by a fac- 
tor of 3 at 100 kpc (see Fig. |3| and there is less apparent 
structure in the ICL component. AGN feedback has been 
quite successful in regulating star formation in the cluster. 
The gas distribution has also been profoundly affected by 
the SMBH model. First, no large, gas rich disc is visible 
in the projected gas density map: the overcooling problem 
have been efhciently removed. We see in Figure [S] that the 
dense unrealistic gas core has disapp eared. When compared 
to the best-fit /3-model proposed bv .Churazov et al] (|2008l ) 
for M87, the agreement is much better. Note that the cooling 
flow, although dramatically reduced, is still present in the 
AGN run, and it can be detected as the density enhance- 
ment within the central 10 kpc of the cluster. Interestingly, 
the knee in the gas density profile seen around 10 kpc in 
our model is also present in the d ata, although much we aker 
and at ~ 5 kpc from the center (|Churazov et al.ll2008l . see 
Fig. 5 in their paper), suggesting the presence of a weak cool- 
ing flow in M87. Another important difference between the 
SF/Quenching runs and the AGN run can be seen at large 
radii in the gas distribution: the gas density in the AGN case 
is 30% larger, showing that gas have been removed from the 
core and stored at large radii (beyond the virial radius) by 
strong shocks similar to the one shown in Figure [2] 
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Figure 6. Effect of the mass resolution on the stellar mass pro- 
files from our AGN model. The blue line i s the stellar mass profile 
of M87 from lGebhardt fc Thomai l l2009t) . 



3.3 The effect of mass resolution 

Using our high resolution simulation with AGN feedback, 
we would like to estimate the effect of numerical resolution 
on our results. We see in Figure [6] the stellar mass profiles 
for the low and the high resolution runs at the final redshift. 
Contrary to what is often claimed in the literature, we do see 
a strong effect of mass (and spatial) resolution in the stellar 
mass distribution. The effect is stronger for the BCG close 
to the center (the stellar mass has increased by a factor of 4 
at 20 kpc) than for the cluster as a whole (the total stellar 
mass within the virial radius has increased only by a factor 
of 2). The effect of mass resolution is smaller for the other 
components (gas and dark matter). 

The strong effect of mass resolution on the star forma- 
tion history of the simulated halo can be interpreted eas- 
ily by comparing the minimum resolved halo mass (opti- 
mistically set to 100 dark matter particles) to the minimum 
mass fo r star forming halos based on atomic cooling argu- 
ments (iGnedinllioOoTlRasera fc Tevssieij l2006l : iHoeft et al.l 
l2006l ). This minimum mass (also referred to as the Filter- 
ing Mass) starts around 10^ Mq before reionization and 
then rises steadily as (1 -I- z)^^^ from redshift 6-7 to the 
final epoch. Resolving this minimum mass before reioniza- 
tion will require a dark matter particle mass below lO'^ Mq, 
a rather strong requirement for cluster-scale cosmological 
simulation. A more flexible criterion based on resolving the 
majority (~ 80%) of star forming halos gives a less stringent 
limit around M^in =i 10® Mq { llie v et al. 2007). Neverthe- 
less, our low resolution run falls short of the corresponding 
required dark matter particle mass by 65, while our high 
resolution run is "only" a factor of 8 above the limit. As ex- 
plained in lLucia fc BlaizotI (|2007l ). BCG are "fundamentally 
hierarchical" objects, that formed their stars very early (80% 
before z — 3) and assembled late (after z — 0.5 in average). 
This effect is directly related to the suppression of cooling 
flows and the associated star formation by AGN feedback. 
This early star formation occurs in rather small mass ha- 
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los l|Lucia fc Blaizotll2007l ) in which star formati on proceeds 
throu gh accretion of diffuse gas in cold streams l|Dekel et alj 
|2009( ). Although BCGs are quite massive objects, it is of 
great importance to resolve properly the earliest epoch of 
star formation , in order to account fo r all the stellar mass in 
these objects. IPuchwein etal\ l|2010l l have reported a simi- 
lar effect in their SPH simulation, although they mentioned 
a significantly smaller effect ('^ 20%) with a different, may 
be more robust feedback model. 

Another interpretation to the rather large resolution 
effect we see in our simulation is an evolution in the effi- 
ciency of AGN feedback. The AGN thermal energy depo- 
sition is performed within a sphere of 4 cells radius. The 
high resolution run will therefore deposit the energy deeper 
into the halo potential well. This might result in a reduced 
overall efficiency . Using the same model than in this paper, 
iBooth fc Schavd liooa) have also reported a rather strong 
dependence of the computed star formation rate on mass 
resolution (see their Fig. 6b). Another solution we would 
like to explore in the future is to recalibrate the AGN feed- 
back model parameters as a function of mass resolution, in 
order to overcome these limitations. 

For both the low and the high resolution simulations, we 
see in Figure |S] that the stellar mass profile has evolved only 
slightly between z = 1 and z — 0. Inside the BCG, we see 
stars expanding slightly, while the stellar halo grows in mass 
more substantially, by almost a factor of 2. This evolution is 
in good qualitati ve agreement with the group-scale simula- 
tion reported bv iFeldmann et ahl (|2010t ). Using our highest 
resolution simulation, we observe that the agreement with 
M87 stellar distribution is very good between 10 and 100 
kpc. It however still deviates quite strongly with observa- 
tions within the inner 10 kpc. We therefore conclude that 
our model with AGN feedback seems to converge from below 
to the correct stellar mass distribution. It is worth stressing 
that the same analysis can be made for our Quenching run, 
and that its converged stellar mass ends up being signifi- 
cantly above the observational limit. From this, we conclude 
that AGN feedback is necessary, not only to regulate star 
formation inside massive galaxies, but also to destroy and 
remove the gas supply in satellite galaxies. 

3.4 The distribution of dark matter 

One important consequence of the overcooling problem is to 
modify significantly the properties of the dark matter halo. 
We have plotted in Figure [7] the projected dark matter den- 
sity for our 3 models, and for the corresponding pure dark 
matter simulation. We immediately see that the dark ha- 
los in the overcooled runs (SF and Quenching) are denser 
and rounder than the pure dark matter case. With AGN 
feedback, we basically recover the halo shape and distribu- 
tion of the pure dark matter case. The effect of gas cool- 
ing on the dark halo has been interpret ed in term of adi- 
abatic contraction of the particle orbits (jBlumenthal et aP 
1 19861 : iGnedin et al1l2004l ) , meaning that individual orbits are 
compressed inward, while conserving the adiabatic invariant 
X = rM( < r). The global shape c hange has been also inter- 
preted bv lDebattista et al.l l|2008h as a transition from boxy 
orbits to more circular ones. These effects have been studied 
quite extensive ly at galactic scales, where they are pro bably 
more relevant l|Abadi et al.ll2009l : IPedrosa et al.]|2009^ . Al- 
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Figure 8. Cumulative total baryonic mass profile (gas -I- stars) 
in the simulated cluster at 2 = 0. The dotted lines are fits to the 
measured profiles for our simple model of adiabatic contraction 
of the dark halo. 
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Figure 9. Cumulative dark matter mass measured in our three 
models at z = 0. Also shown in blue is the profile we obtained in 
the dark matter only (DMO) simulation. In each case, the dotted 
line corresponds to our analytical model of adiabatic contraction. 
Note that in the AGN feedback case, we see a slight adiabatic 
expansion of the dark halo. 

though we are dealing with a much larger object, we recover 
very similar properties for our dark halo, because of over- 
cooling. We have plotted in Figure [9] the cumulative dark 
matter mass profiles for our three runs, plus the Dark Matter 
Only (DMO) simulation. The DMO profile has been multi- 
plied by 85% to allow a direct comparison with the baryonic 
runs: it is fitted with an accuracy better than 5% down to 
10 kpc by a NFW profile with a concentration parameter 
c — 7.5. The fit is shown as the blue dotted line on the same 
figure. The dark matter profiles for the SF and Quenching 
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Figure 7. Maps of the dark matter overdensity distribution in tlie simulated cluster at 2 = for our 3 models. Also shown for comparison 
is the dark matter density map we have obtained in a dark matter only (DMO) simulation of the same cluster. 



runs are very similar, except in the very centre of the cluster. 
They all appear significantly adiabatically contracted. Inter- 
estingly enough, the AGN case appears slightly expanded, 
when compared to the DMO simulation. We will now use 
AC theory to explain these trends. 

Gnedin et al. hl2004D have revisited the original paper 



of iBlumenthal et aL ( 19861 ) on AC theory, stressing that 



the original assumption of purely circular orbits was leading 
to an overestimate of the baryon-induced dark halo con- 
traction. They presented a numerical implementation for 
their modified AC theory in which the particle orbit distri- 
bution was allowed some radial components and predicted 
the contracted dark matter distribution, given the initial 
dark matter profile and the final baryonic distribution. We 



present in the Appendix a simple analytical model to ac- 
co unt for the ad i abatic contraction of the dark halo, based 
on lCnedin et all ^2004 ) theory. 

In Figure [8] we have plotted the total baryonic mass 
Mbar as a function of the final radius for our 3 different mod- 
els. Although the actual distributions are different from our 
simple model, we have fitted them using a constant surface 
density, truncated disc of mass md and size r^. The fits are 
very similar in the SF and Quenching case, so we have used 
only one model with ~ 5 x 10^^ M© and ~ 20 kpc. 
These values, consistent with the mass of our simulated BCG 
(including the concentrated gas component in the Quench- 
ing case) are much larger than any observed BCG in ha- 
los of similar masses, and are again a manifestation of the 
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overcooling problem. In this case of strong baryonic concen- 
tration, our analytical model with <C Va applies (see the 
Appendix): we have plotted the corresponding AC theory 
prediction in Figure [9] as th e dott ed line, showing convinc- 
ingly that the lGnedin et all l|2004) model, in our simplified 
formulation, works quite well: the fit is better than 20% 
down to 10 kpc in radius. Within 10 kpc, our fit is not as 
good. We believe that since we are entering the scale of both 
the disc and the dark matter inner region (r < rs), the adi- 
abatic contraction model doesn't apply anymore. When the 
halo first forms and virializes at high redshift, the baryon 
fraction in star and gas during violent relaxation is likely 
to play an important role, explaining why the dark matter 
profiles in the SF and Quenching runs differ below 10 kpc, 
although the total baryon mass profiles are very similar at 
redshift zero. 

The AGN feedback model gives us a much more ex- 
tended baryonic mass profile. We have fitted it using the 
constant surface density truncated disc model with param- 
eters nid ~ 2 X 10^^ M0 and rd ^ 700 kpc. We can see 
in Figure [8] that our fit is far from being perfect (especially 
within the BCG), but it captures roughly the total baryon 
distribution, which is now mostly dominated by the hot ex- 
tended gaseous halo. The values for and ra are now closer 
to the total mass and total size of the hot halo than those of 
the BCG. In this case, the AC prediction cannot be worked 
out analytically. We have therefore solved for the real root 
of the third order polynomial defined by Equation lAll and 
plotted the r e sult i n Figure [9l Again, the prediction from 
ICnedin et al.l (|2004l ') theory is very close (within 20% above 
10 kpc) to the measured dark matter profile. In the AGN 
feedback case, we see that the dark matter has expanded 
slightly, when compared to the pure dark matter case, and 
that this "adiabatic expansion" appears to be well captured 
by the same adiabatic invariant for the orbits of dark mat- 
ter particle. Note that this expansion is very small, as it 
affects the dark matter mass distribution by less than a few 
percent. We have also checked that us ing the complete nu- 
merical solution of ICnedin et al.l (|2004l ) in conjunction with 
the simulated baryon mass profile (instead of our simple disc 
model) does not affect our conclusions. 

A good diagnostic of how much concentrated the simu- 
lated halo should be to match the observations is to compare 
the total mass profile with the one derived from dynamical 
arguments using M87 optic al and X-ray data (see Fig. llOp . 
ICebhardt fc Thomas! ()2009l l have derived a complex mass 
model for M87, including the presence of a central SMBH, 
stellar BCG and dark matter halo. Their mass model is 
compared to our total mass profiles from our three differ- 
ent models in Figure [TOl The SF and Quenching models are 
too much concentrated and overestimate the total mass by 
a factor of 2 to 4 in the inner 100 kpc. The AGN model is 
in much better agreement, although slightly below, for both 
the low and high resolution runs. 

Another diagnostic on the total mass distribution is the 
temperature profile of the hot. X-ray emitting gas. We have 
plotted in Figure [TT] the mass- weighted temperature profile 
for our various simulations. The temperature at 100 kpc is 
3 keV in the AGN cas e, in good agreement w ith the observed 
value of 2.8±0.2 keV (jChurazov et al.ll2008l l. Without AGN 
feedback, we obtain a larger temperature (around 6 keV) be- 
cause of the higher mass concentration. If one looks at the 
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Figure 10. Cumulated total mass profile (baryons + dark mat- 
ter) in the simulated cluster at 2 = 0. The blue line in the is 
the mass profile deduced from stellar kinematics and X-ray data 
IGcbhardt fc Thomas 200a) . 
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Figure 11. Mass- weighted gas temperature profile in the sim- 
ulated cluster at 2 = 0. The blue line in the is the temperatur e 
profile inferred from X-ray data in Virgo llChurazov et al.1l2008l ). 



temperature profile, we do not get a good match with X-ray 
data, even including AGN feedback. We are slightly hotter 
than the observations between 10 and 100 kpc, and within 
10 kpc, we see a sharp drop of temperature, due to cool- 
ing gas feeding the central galaxy and sinking towards the 
central AGN. We have checked that hydrostatic equilibrium 
is satisfied down to 10 kpc. Matching X-ray temperature 
profiles will probably requires additional physics, such as 
cosmic rays propagation and magnetic fields, and/or a more 
realistic AGN feedback scheme. 
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Figure 12. Cumulative (top panel) and local (bottom panel) 
baryon fractions at redshift 2 = in the simulated cluster. The 
universal baryon fraction is shown as the blue horizontal line, 
while the virial radius is indicated by the dashed vertical line. 
Each model is labelled using the standard color scheme. In the 
AGN feedback case, we observe a slight deficit (~10%) of baryon 
within the virial radius. Using the local baryon fraction, we see 
that these missing baryons are located in between R^jir and 2i?^i,. . 
In the upper panel, the cumulated gas fraction is also shown as a 
dashed line. 



3.5 The baryon fraction 

The most important consequence of the disappearance of 
the overcooling problem is that we expect the baryonic mass 
distribution to be in much better agreement with observa- 
tional constraints. Using X-ray data, it is indeed possible to 
estimate the gas profile in large clusters, while the stellar 
mass can be measured using optical data. iGonzalez at al.l 
l|2007t ) have computed the baryon mass fraction within rsoo 
for a large sample of groups and clusters. They have found a 
slight deficit of baryons /bar = 0.133±0.004, when compared 
to the universal baryon fraction as measured by WMAP 
Qb/ilm = 0.176 ± 0.008 (|Spergel et al.ll2003l . bOOTl ). From 



Figure 1121 we see that the SF and Quenching runs show a 
slight baryon excess with /bar — 0.16 while our universal 
baryon fraction ^It/i^m was set to 0.15 in our simulation 
(see Table [2l). This is a direct consequence of overcooling 
(jKravtsov et al.l I2OO5I ). On the contrary, the AGN model 
show a clear baryon deficit with /bar — 0.13, in striking 
agreement with the observed average value. Note however 
that we used a universal baryon fraction lower than the 
WMAP estimate, so that we cannot claim that our model is 
fully consistent with the data yet. Nevertheless, this baryon 
deficit in our model is the consequence of the repeated effect 
of AGN-driven shocks and convective motions, pushing gas 
outside the virial radius. One can see from the local baryon 
fraction profile in Figure [12] that this gas accumulates in a 
region between 1 to 2 virial radii around the cluster, be- 
yond which the cumulated baryon fraction converge to the 
universal one. 

Our SF model compares very well to the AMR simula- 
tion performed by iKravtsov et al] l|2005l l. showing a slight 
excess of baryon at the virial radius. From these standard 
galaxy formation models, we obt ain gas properties that com- 
pare favorably to X-ray data (|Nagai et al.l |2007| ). We see 
indeed in Figure [12] that in the SF model, the gas mass 
fraction is slowly decreasing towards the center, with a sig- 
nificant deficit at the virial radius. This behavior is tradi- 
tionally explained by the join t effect of cooling and star for- 
mation (|Voit fc Brvanll200ll ). The price to pay is however 
to form too many stars and cold gas i n the cluster, as con- 
firmed by previous numerical models (iBorgani et al.l l2004l : 
iKravtsov et al.ll2005l : iBorgani fc Kravt sov 2009). Our sim- 
ple Quenching model is making things better for the stellar 
component, but now the gas mass distribution is too concen- 
trated (see Fig |12|) . Only with AGN feedback can we obtain 
a small stellar mass fraction and in the same time a small 
gas fraction. We note however that, in our AGN model, the 
gas mass profile is not as steep as suggested by X-ray data. 
We could probably reproduce the gas profile of the SF model 
using another, more efficient AGN feedback scheme. 

Recently, jPuchwein et alj (|2010 l) have simulated a large 
number of groups and clusters with the SPH code GAD- 
GET, using a mass resolution that is only about a factor 2 
lower than ours and a spatial resolution of 2.5kpc (compared 
to 1 kpc at low res. or 0.5 kpc at high res. here). Nevertheless, 
they also found that with AGN feedback, the total baryon 
fraction was below the universal value. More interestingly, 
using a high universal baryon fraction {Qt/^lm — 0.165), 
they report a stellar mass fraction of /* ~ 0.05, quite inde- 
pendent of the parent halo mass. In our case, for our Virgo- 
like cluster with A/vir — lO^** Mq, we obtained /, ~ 0.01 in 
the low resolution case and /* ~ 0.02 in our high resolution 
simulation. Our value is in bet ter agreement w i th the obser- 
vational estimate proposed byl Lin e t al. ' ('2003'. '2004). while 
the higher value found by Pu chwein et a l. (2 010.) is in better 
agreem ent with the observations repor ted in Gonzalez et al.l 
(j2007l ). Using the COSMOS survey, iGiodini et ahl (|2009D 
have estimated stellar mass fractions for a large sample of 
galaxy clusters and groups. For our simulated halo mass, 
M500C ~ 8 X 10^'' M0 , they report stellar mass fraction rang- 
ing from 2% to 5%. U sing the original feedba ck model of 
iBooth fc Schavd (|2009D in the GADGET code, iDuffv et all 
(20l3) have also computed the predicted baryon and stellar 
mass fraction of a large sample of groups extracted from 
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a cosmological simulation. Although they report a simi- 
lar baryon deficit within the virial radius, they obtained 
/t ~ 0.03 for an even larger universal baryon fraction 
= 0.18. We see that there is a consensus about a 
strong reduction of the stellar mass fraction in groups and 
clusters thanks to AGN feedback. The extent of this reduc- 
tion seems to depend quite sensitively of the details of each 
implementation, and possibly to the nature of the code (SPH 
versus AMR). As suggested by observations also, we note 
that the exact stellar and baryon fraction probably varies 
from halo to halo. We would like also to stress that all the 
reported simulations, including ours, are probably not fully 
converged yet. 



4 SUMMARY AND CONCLUSIONS 

We have simulated the formation of a Virgo-sized galaxy 
cluster to study the effects of feedback on the overcooling 
problem. The impact of AGN feedback on the distribution 
of the baryonic mass is strong, and in good agreement with 
previous SPH simulations: star formation in massive galax- 
ies is drastically reduced. At the same time, left-over gas 
is very efficiently removed from the core of the parent ha- 
los, where it would have otherwise accumulated. In order to 
quantify the effect of AGN feedback, we have run two other 
reference simulations: one model with only star formation 
and supernovae feedback (the standard scenario) and one 
model for which we have artificially prevented star forma- 
tion to occur in massive enough spheroids (the quenching 
scenario). 

A detailed comparison of the three models clearly 
demonstrates that AGN feedback is needed to control star 
formation in the central BCG, but also to unbind the over- 
cooling gas from the cluster core. We also clearly identify 
the effect of the baryon dynamics on the dark matter mass 
distribution on large scale. Interestingly enough, in case 
of AGN feedback, we observe the adiabatic expansion of 
th e dark halo, an effe ct well modeled by the AC theory 
of iGnedin etld] (|2004 ). A comparison of our simulation re- 
sults with observational data for Virgo and its central galaxy 
M87 rules out the standard model, but also the quenching 
model. On the contrary, our simulation with AGN feedback, 
although not fully converged yet, shows a much better agree- 
ment with M87 data in term of mass distribution. In particu- 
lar, we obtain a significantly reduced baryon fraction within 
th e virial rad ius, in agr eement with observatio ns compiled 
bv lLin et all (2003) and iGonzalez et al.l ^200l\ ). We clearly 
identify in our simulation that gas is removed from the core 
of the cluster by convective motions and/or strong shocks, 
and accumulates in a region just outside the virial radius. 
When compared to gas profiles inferred from X-ray data, 
our AGN model produces too shallow gas distribution, sug- 
gesting that we probably need even more powerful feedback 
processes. 

Our cluster formation simulations with AGN feedback 
have not fully converged yet - as we increase the resolution, 
we find a stellar mass profile for the BCG that is in bet- 
ter agreement with observations, but it is still too low by 
about a factor of two. We are still missing the lowest mass 
galaxy population, which could provide the missing stellar 
mass in the central elliptical galaxy. We also note that, in the 



current picture, AGN feedback is a morphologically depen- 
dent process: it only directly affects galaxies with SMBHs, 
i.e. galaxies with a significant bulge/spheroid component. 
Higher resolution studies would be needed, in order to reli- 
ably model this distinction, so that star formation in disky 
galaxies is not artificially suppressed. 
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APPENDIX A: ADIABATIC CONTRACTION 
MODEL 

If one defines the initial radius of each dark mat ter shell as Vj 
and its final, adiabatically contr acted value rf lAbadi et alj 
(i2009D have proposed to capture ICnedin et al.l (|2004 ) model 
using the following simplified model 



1 + a 



Mf 



with 



0.68. 



(Al) 



The original fBlumenthal et ahl (Il986l ) can be recovered using 
a — 1. The final cumulated mass distribution is computed 
using 

Mf = Mi^{rf) + Kh^,{rf) = f^^M^in) + M^^,{ri) (A2) 

where we have assumed that the initial dark matter mass is 
conserved during AG. The initial dark matter distribution 
is described using the analytical NEW profile 



M^n) = M200 



log(l + x) -x/il + x) 
log(l + c) - c/(l + c) 



(A3) 
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where x = Vi/rs and = r2oo/c- M200 is the total vhial 
mass. For the baryonic distribution, we assume a constant 
surface density disc with size rd and mass md, so that 



The dark matter mass fraction is computed using fd = 
1 — md/M2oo- The model we considered in Equation I A4I for 
the baryonic mass distribution has been chosen that simple 
on purpose: inserting Equation IA4I into the AC relation in 
Equation lAll one clearly sees that we have to find the only 
real root of a third order polynomial equation with unknown 
Tf/ri. This can be done quite easily with any root finder. In 
case the disc size is small enough (namely if rd <g r^), the 
AC model is fully tractable analytically by noticing that for 
a; ^ 1, one has Mi oc x^. We therefore have 



where the constant can be determined by continuity. 




(A4) 
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— ~ constant for rt < rd 
ri 
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